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While adaptive optical systems are able to remove moderate wavefront 
distortions in scintillated optical beams, phase singularities that appear in 
strongly scintillated beams can severely degrade the performance of such an 
adaptive optical system. Therefore, the detection of these phase singularities is 
an important aspect of strong scintillation adaptive optics. We investigate the 
detection of phase singularities with the aid of a Shack-Hartmann wavefront 
sensor and show that, in spite of some systematical deficiencies inherent 
to the Shack-Hartmann wavefront sensor, it can be used for the reliable 
detection of phase singularities, irrespective of their morphologies. We provide 
full analytical results, together with numerical simulations of the detection 
process, (c) 2008 Optical Society of America 
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1. Introduction 

The Shack-Hartmann wavefront sensor (SHWS) is a widely used system to measure the 
shape of the wavefront of an optical beam that was scintillated after propagating through a 
turbulent atmosphere. The reasonable simplicity of the SHWS makes it ideal for use in an 
adaptive optical system, which is to correct the wavefront distortions in such a scintillated 
beam. 

Various techniques exist^"^ to compute the shape of the wavefront from the raw data 
received as output from a SHWS. However, in strongly scintillated beams these techniques 
tend to faiF due to the presence of phase singularities, which represent optical vortices^' ^ in 
the propagating scintillated beam. The failure of these techniques is to some extent a result 
of the apparent insensitivity of the SHWS to the phase function of a phase singularity, 
which led to it being referred to as hidden phase. Phase singularities are points where the 
phase is undefined and where the phase around the singularity goes continuously through all 
phase values from to 27r. The direction (right-handed or left-handed) of increasing value 
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indicates the topological charge (±1) of the phase singularity. The phase function around 
a singularity can in general be anisotropically scaled along a particular direction, which 
gives rise to its morphology. An isotropic singularity is referred to as being canonical and 
anisotropic singularities are called noncanonical. 

Initial least-squares phase reconstruction methods/"^ that were used in adaptive optical 
sytems, were based on a model for the phase slopes that makes the assumption that the 
phase function is continuous. The existence of phase singularities violates this assumption. 
Therefore, these methods can only be applied in application where the phase distortions are 
not large enough to generate phase singularities. 

More recently, various authors^^"^'^ considered the separation of the phase gradient function 
into a continuous (curl- free) part and a solenoidal or rotational part. The former can be 
used in the adaptive optical system to correct the continuous part of the phase distortions. 
The latter is generally treated as a noise term, which is discarded. Since it represents the 
contribution of the phase singularities, the solenoidal part of the phase distortions cannot 
be removed in the same way that the continuous part of the phase distortions is removed. 
Neither does it simply go away after the continuous part has been corrected. It is therefore 
important to consider this solenoidal part of the phase distortions more carefully in order 
to find a way to get rid of the phase singularities. In this paper we consider the detection of 
the locations and topological charges of the phase singularities with the aid of the solenoidal 
part of the phase gradient function. We show that the SHWS is not completely insensitive 
to the phase of a singularity. In other words, this phase is not completely hidden, but is 
contained in the solenoidal part of the output of the SHWS. One can therefore use it to 
locate the phase singularities that are present in the wavefront. 

The detection of the phase singularities is significantly influenced by the averaging effects 
of the SHWS subapertures. This averaging effect has not to date received much attention. 
Herrmann^^ mentioned the model error due to the averaging gradients, but did not give a 
detailed analysis. Aksenov^^ calculated the average wavefront slopes with a regularization 
procedure that produces continuous infinitely differentiable functions, which is not true when 
there are singularities in the phase function. In this paper, we analyze the effect of the 
averaging process on a phase singularity, by performing the actual integration of the phase 
gradient function of such a singularity over the area of a subaperture. From this result one 
can then show that it is possible to extract the required information of the phase singularities 
regardless of the effects of the averaging process. For a continuous vector field the singularities 
can be identified by computing the curl of this vector field. The output from a SHWS is in 
the form of a sampled vector field that approximately represents the gradient of the phase 
function of the scintillated beam. We implement the curl operation that operates on this 
sampled vector field in terms of finite differences. This finite difference curl operation is 
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referred to as the circulation of the sampled vector field. The complete detection process 
thus consists of the averaging process performed by the SHWS and the circulation process 
performed on its output. The effect of the morphology and relative position of the singularity 
on this detection process is also investigated. 

The remainder of the paper is organized as follows. In Section [2] we'll discuss the basic 
principle upon which the detection of phase singularities is based and provide some notation. 
The operation of the SHWS is discussed in Section [31 In Section H] we discuss the detection 
of phase singularities with the output of the SHWS. We analyse this for a canonical phase 
singularity and compute the peak value in Subsection 14. Al and then repeat this analysis for 
the more general noncanonical phase singularity in Subsection 14. B[ The complete expression 
of the circulation for a noncanonical phase singularity is given in Appendix El A numerical 
simulation is provided in Section [5] to test the procedure on more realistic data. A summary 
and conclusions are provided in Section [61 

2. Basic principle of operation 

The phase function of a scintillated optical beam can in general be represented by the sum 
of a continuous phase function and an arbitrary number of phase singularities, 

e{x, y) = edx, y)+J2 -Xn,y- yn] Q!n, /^n), (1) 

n 

where 6c{x, y) is the continuous phase function and (j){x — x„, y — yn', an, Pn) represents the 
phase function of a noncanonical singularity, located at {xn,yn)- The latter phase function 
is given by, 

^(x + iy) + C{x-iy) 



1 

{x,y]a,f3) = --In 



^* {x — iy) + (* {x + iy) ' ^ ^ 

where * represents the complex conjugate and the morphology of the singularity is parame- 
terized by, 

^ = cos(a/2)exp(z/3/2), 

C = sin(a/2)exp(-^/5/2), (3) 

in terms of the morphology angles < a < vr and < /3 < 2tx. 

The presence of a singularities in a phase function can be determined with the aid of a 
closed line integral over the gradient of the phase function, 

Ve-dl = r27r, (4) 

c 

where C denotes the closed integration contour; V9 is the gradient of the phase function 
and r is an integer that presents the net topological charge of all the singularities enclosed 
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by the contour. Unfortunately the closed line integral in Eq. (jl]) is not convenient for the 
practical implementation of singularity detection. 

To implement a singularity detection procedure with a physical system one can use the 
output of a SHWS, which measures the gradient of the phase function. For the moment we'll 
assume that this measurement is ideal and we'll return to the effect of the finite subapertures 
of the SHWS in Section [31 One can view the output of a SHWS as a vector field G = 
Gx{x,y)x + Gy{x,y)y. If this vector field represents the gradient of a continuous phase 
function, we must have, 

dG^{x,y) ^ dGy{x,y) 
dy dx 

If one finds that this is not the case then it indicates that G cannot be the gradient of a 
continuous phase function. This points to the presence of phase singularities. The extent 
to which the equality in Eq. ([5]) fails is given by the difference between its left-hand and 
right-hand sides, which becomes the curl of the vector field. Since the vector field only has 
X- and y-components that only depend on x and y, its curl is a scalar which represents the 
2;-component given by 

D(x,y)=VTxG(x,y). (6) 

If the phase function of a scintillated optical beam is continuous (without phase singular- 
ities), the curl of the gradient of the phase function would give zero. On the other hand, if 
there are phase singularities in the phase function of the beam then the curl of this gradi- 
ent is not zero.^^ Each phase singularity gives a Dirac delta function at the location of the 
singularity, 

VxV0(x, y, a, P) = T2Tx5{x)5{y), (7) 

where r(= ±1) is the topological charge of the phase singularity. So the curl of the gradient of 
a phase function 6{x, y) that contains phase singularities gives a sum of Dirac delta functions, 
each multiplied with the topological charge of the singularities, 

VxV^(x, y) = 2^Y. ^n5{x - Xn)S{y - (8) 

n 

By the same token the curl of the output vector field of a SHWS, shown in Eq. ([6]), gives 
us, according to Eq. ([8]), a distribution of singularities, showing us where they are located 
and what their topological charges are. One can easily compute D{x, y), as shown in Eq. ([6]), 
from G and then use the result to identify the phase singularities in the wavefront. 

In the above discussion we assumed the ideal case, which ignores the effects of the discrete- 
ness and finite size of the subapertures of the SHWS. In the following section we consider 
the SHWS more carefully. We'll see that the subaperture size sets a scale beneath which 
singularities become unobservable and it also modifies the weight of 27r with which each 
singularities is multiplied. 
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3. Shack-Hartmann wavefront sensor 

The SHWS uses a lenslet array to sample the wavefront, as shown in Fig. [H The slope of 
the wavefront for each of these samples is given by the location of the focal point formed 
by each lenslet in its back focal plane. To simplify our analysis we'll assume that the lenslet 
array is a square array of lenslets, each with a square shaped subaperture, as shown in Fig. [2j 
A detector array is placed in the back focal plane of the lenslet array. A small subarray of 
this detector array is dedicated to each lenslet of the lenslet array. This subarray is used to 
determine the location of the focal point produced by each lenslet. The location of the focal 
point for the (m, n)-th subaperture (lenslet) is given by a position vector u™'" that is given 
by the centroid of the measured intensity distribution in the back focal plane, 

„m,n /ff -^(u)u d^M ^ 

where uq"^'" would be the location of the focal point for a normally incident plane wave; 
/(u) is the intensity distribution over the detector plane; H is the window of the subarray 
of detectors on the detector plane; and d'^u is the two-dimensional integration measure on 
the detector plane. 

The average phase slope for each subaperture is given by the location of the focal point, 
as determined by the intensity centroid in Eq. This relationship follows from the Fourier 
relationship, as portrayed in Fig. [3], between the phase tilt of the incident wave in front of a 
lens and the location of the resulting focal point behind it. So the average phase slope over 
the (m, n)-th subaperture can be expressed as 

rim.n _ /nV^(x) d^x _ k 

- In - 7" ' ^''^ 

where V^(x) is the gradient of the phase function of the incident wave; x denotes the two- 
dimensional position vector on the lenslet plane; k{= 2tt/X) is the wave number; / is the 
focal length of the lenslets; Q is the area of the lenslet subaperture; and u™'" is the location 
of the (m, n)-th focal point on the detector plane, as defined in Eq. ([9]). The G™'""- values 
given by Eq. (ITUl) represent the sampled output vector field of the SHWS. Each value is 
associated with a point in the center of the particular subaperture, as denoted by the dots 
in Fig. H 

4. Detection of phase singularities 

In the previous section we saw that the output of the SHWS in a practical adaptive optical 
system is a sampled vector field G with each sample representing the averaged phase slopes 
over one subaperture. The solenoidal part of this vector field contains the information about 
the locations of the singularities, as pointed out in Section[2l To extract this information from 
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the sampled vector field, one needs to implement the curl operation of Eq. ([6]) numerically. 
Such a numerical implementation is equivalent to the line integral of Eq. (jlj). The numerical 
computation is done by computing what we refer to as the circulation, given by. 



where w is the subaperture window size and the superscript m, n denotes the subaperture 
index. The physical implementation of the circulation process is presented in Fig. |H The 
samples used in the calculation are denoted by the points in the centers of the subapertures 
shown in Fig. |H The circulation represents a line integral performed over the four subaper- 
tures along a contour denoted by the dashed-line. Note that the result of this circulation 
operation represents a value D that should be associated with a point in the center of the 
four subapertures. However, we are assigning it the same index m, n that is associated with 
the upper left subaperture. 

The result of this computation ]J^'^ is a distribution of the topological charges of the phase 
singularities. It is positive (negative) at the locations of phase singularities with positive 
(negative) topological charges and it should be zero where there are no phase singularities. 
However, this distribution is affected by the averaging process, the sampled nature of the 
data and, of course, noise. As a result the values are not exactly zero when there are no 
singularities. Moreover, the value of Z}™'" at the location of a single positive singularity is 
not 271 as one would expect if the circulation is an exact implementation of the line integral 
Eq. (j4]). In the next two subsections, we analyze the effects of the practical implementation 
of the singularity detection procedure. 

4. A. Canonical phase singularities 

Here we consider the situation where a canonical singularity is located in the center of the 
four subapertures, as shown by point A in Fig. HI We define the origin (0, 0) of our coordinate 
system at this point. The complex amplitude function of a canonical phase singularity at the 
origin can be expressed as {x ± iy)/r = exp(±i0), where the sign indicates the topological 
charge of the singularity and r and are respectively the radial coordinate and the azimuthal 
coordinate. Here we consider a positively charged singularity. The phase function of the 
singularity is simply the azimuthal coordinate (p, provided that the singularity is located at 
the origin. The gradient of the phase function of a canonical phase singularity can then be 
expressed in Cartesian coordinates, as. 



D 



,m,n 





V(f){x,y) 



xy — yx 
x"^ + 



(12) 
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Substituting Eq. (fT2!) into Eq. (fTOj) . we obtain an analytical expression for the average phase 
slope in the (m, n)-th subaperture, given by, 







J-w 


J —W 




In 2 


\4:W 


2^ 



i-y). (13) 

The integration boundaries are determined by the location of the (m, n)-th subaperture 
within the four subaperture area shown in Fig. HI Due to the symmetry of the phase func- 
tion of a canonical singularity at the origin, the average phase slope values for the other 
subapertures will be the same apart from a possible change in sign. Then, according to 
Eqs. (mD and ([T3D, the value for D™'" will be, 

D'"'" = 4m;G'^'" = 71 + 2 ln(2) = 4.527887. (14) 

We note that the value of the circulation is vr + 2 In 2 and not the 2tt that one finds for 
the analytical case given by Eq. (jl]). It is the averaging process that is responsible for this 
difference and not the finite differences of the calculation process. One finds that if the same 
finite difference circulation calculation, Eq. (fTTj) . is performed on a sampled gradient function 
of a canonical phase singularity without the averaging process the result is indeed 2-7r. 

Although it is the averaging process that causes the difference, we note that this circulation 
value is independent of the subaperture window size w. This is because the phase function of 
a singularity is scale invariant. So the difference is produced simply because some averaging 
takes place, but it does not matter how big the subapertures are over which the slopes are 
averaged. 

When the singularity is not located at the origin, in the center of the four apertures, but 
at some other location such as the position B at (xq, yo) in Fig. HI we can expect the value of 
the circulation to change. The expression for the gradient of the phase function now becomes 

v7Af \ (x - xo)y -{y- yo)x 

V^(.-.„,V~V„) = j—^,^^—^. (16) 



Using Eq. f[T5|) and Eq. fjTOj) . one can compute the average phase slopes G for the four 
subapertures. Then, with the aid of Eq. f|TT]) . one can compute D"^'^. The final expression 
is rather complicated. It can be obtained from the more general expression for arbitrary 
morphology provided in Appendix [XJ as explained there. 

The circulation D is shown in Fig. [5] as a function of the location of the singularity in 
terms of normalized coordinates fi = xq/w and u = yo/w. A topview of this function over the 
region —2 < /i < 2 and —2 < z/ < 2 is shown in Fig.[5](a). The precise shape of the circulation 
function as one- dimensional functions of /j, are shown in Fig. [5](b) for u = (middle line) 
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and for u = fi (diagonal). At the origin (/i, i^) = (0,0) the circulation function has a peak 
with the value D = 4.53, consistent with Eq. dHj). Away from the origin the value of D 
decreases rapidly and approaches zero for > 1. In some regions the value of D drops 

below zero. At the corners of the four subaperture area, where |z/| = 1, the circulation 
function forms negatively valued peaks with the value, 

= I - 2 ln(5) - 2 arctan(2) + 5 ln(2) = -0.396641. (16) 

The fact that the circulation function never becomes exactly zero implies that the existence 
of a singularity at some point in the output plane gives nonzero circulation values in all other 
parts of the output plane. 

In view of the fact that the peak value of the circulation function is vr + 2 ln(2) instead of 
27r, it is interesting to note that, if one would add the four sample values of Z}"*'" closest to 
the location of a singularity - i.e. the four values that surround the singularity - then the 
result is closer to 27r. 

4.B. Noncanonical phase singularities 

Singularities in strongly scintillated beams are in general noncanonical. It is therefore nec- 
essary to know how the computation of the circulation is affected by the morphology of the 
singularity. The phase function of a singularity with an arbitrary morphology is given by 
Eqs. and ([3]). The gradient of this phase function can be expressed by, 

= x^{l + A)-2yxB + y^{l-Ay ^'^^ 

where, 

A = sin(a)cos(/3), (18) 

B = sin(a)sin(/3), (19) 

C = cos(a). (20) 

with the morphology angles a and (3, as defined in Eq. ([3]). The singularity can be translated 
to any location (xq, yo) by replacing x ^ x — xo,y ^ y — yo in Eq. (fTTl) . 

We use the same procedure that was use in Subsection 14. Al to analyze the noncanonical 
case. We use Eq. (flOl) to compute the output of the SHWS for the phase gradient given 
in Eq. f[T7|l . shifted to (xo,yo)- Then we substitute the result into Eq. flTT]) to compute the 
circulation D. In Appendix |A] we provide the complete expression of this circulation function 
for a singularity with an arbitrary morphology, located at an arbitrary point (in normalized 
coordinate n = Xq/w and u = y^/w) within the four subapertures, shown in Fig. HI 

In Fig. [6] we show the peak value at (/i, v) = (0, 0) for D as a function of the morphological 
angles a and (3. Note that the peak is positive (negative) forO < a < 7r/2 (7r/2 < a < 
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it). At a = 7r/2 there is a discrete jump, which represents the change in the topological 
charge of the singularity. The fact that there is such a large difference between the values 
on either side of the jump, indicates that the circulation D, calculated from the output 
of SHWS, can in principle determine the topological charge of a singularity even if it is 
severely anisotropic. However, in a practical situation, singularities tend to have such severe 
anisotropic morphologies only when they appear in oppositely charged pairs close to each 
other. In such situations their circulation values will in general partially cancel each other, 
making them difficult to identify. Near the jump at a = 7r/2 the value of D has its greatest 
fluctuation as a function of (3. Next to the jump D fluctuates between vr and 2-7r. When a 
approaches the canonical values of or vr the value of D tend towards its canonical peak 
value of 4.53 and the fluctuations as a function of /5 diminish. 

In Figs. [7] and [8] we provide plots of D over the region — 2 < /i < 2 and — 2 < z/ < 2, for 
two different morphologies of the singularity. First we consider the case when a = tt/A and 
(3 = 71. This represents a singularity with a moderate anisotropy oriented along the y-axis. 
The topview of the circulation function for this case is shown in Fig. Wi^) and the shape of 
the circulation function is shown in Fig. H^b) in terms of three one-dimensional functions: 
along 1/ = (/i-line), along fi = (z/-line) and along the line where u = fi (diagonal line). 
The peak value at the origin is about 4.5. The function then decreases toward zero away 
from the origin. We note that the shape of the circulation function is more anisotropic than 
the shape in Fig. [5] in that the respective widths of the peak along the fi- and t^-directions 
are not equal. There are still regions where the function becomes negative. 

Next we consider the case where a = An/d and (3 = 7r/2. This represents a highly 
anisotropic singularity oriented diagonally along the line where u = —fi. We show the topview 
of the circulation function for this case in Fig. [7](a) and the shape of the circulation function 
in Fig. [71(b) in terms of three one-dimensional functions: along u = (middle line), along 
the line where u = fi (diagonal-I line) and along the line where u = —fi (diagonal-II line). 
The peak value at the origin is now about 5. The function decreases away from the origin, 
but the decrease is much slower along the orientation of the singularity. The shape of the 
circulation function is therefore much more anisotropic. 

5. Numerical simulation 

Here we present a numerical simulation to test the singularity detection procedure that was 
analytically investigated in the previous sections. We simulate the propagation of a Gaussian 
beam over a distance of 100 km through a turbulent atmosphere. We use the well known 
numerical method of Ref. 17 to perform the simulation. The strength of the turbulence is 
parameterized with = 4 x 10~^^ m^/"^ and we use 10 equally spaced phase screens to 
simulate the 10 km thick turbulent layers. This method provides a reasonable agreement 
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between real world data and the simulation data. The phase of the beam is distorted when 
it reaches the system aperture and in our example we find four phase singularities in the 
wavefront, as shown in Fig. ^a). This beam then passes through our simulated SHWS, 
from which we obtain the sampled vector field G, computed with Eq. ffTOj) . The circulation 
function D is then computed with Eq. f|TT|) . The resulting circulation function is shown in 
Fig. MJo). Two of the singularities, one positive and one negative, are easily identified from 
their respective circulation values of 2.95 and —3.54, in the lower-left corner of Fig. M^a). 
Note that the magnitudes of both these values are smaller than 27r. This would be due to 
a combination of the fact that the singularities have noncanonical morphologies, the fact 
they are not located at the ideal location in the center of four subapertures and noise that is 
present in the phase function. Integrating over a 3 x 3 neighborhood around these singularities 
we obtain values of 6.20 and —5.86, respectively, which are closer to ±2tt. 

The other two singularities in the upper-right corner of Fig. M^a) are much closer to 
each other. Therefore, their individual circulation functions overlap and, having opposite 
topological charges, they partially cancel each other. As a result the circulation peaks that 
represent these singularities are severely diminished. For example, the positive peak for this 
pair of singularities in Fig. Mjo) has a value of only 1.41. Oppositely charged singularities 
that are located closer to each other are therefore more difficult to detect. 

6. Conclusion 

The phase gradient that is produced as output of a SHWS contains information about the 
continuous phase function of the incident wave, but also of the phase singularities in the 
wavefront. A least-squares projector^^"^^ can be used to extract the information about the 
continuous phase and can be used to correct contious phase distortions. The information 
about the singularities is contained in the solenoidal part of the phase gradient. The curl 
of this part gives a topological charges distribution, which represents the locations and 
topological charges of the singularities. Theoretically each positive (negative) singularity 
should be indicated by a value of 27r (— 27r) in the topological charge distribution. 

The averaging process inherent to the SHWS has a significant effect on the computed 
topological charge distribution. Instead of the theoretical value of 27r, the actual value that 
is produced is at most about 4.53. The precise location of the singularity relative to the 
subapertures in the SHWS, as well as the morphology of the singularity produce further 
variations in the value of the topological charge distribution at the location of a singularity. 
Nevertheless, these values are generally large enough to identify an isolated singularity. It 
is therefore possible to extract the information of the location and topological charge of the 
singularities from the output obtained from a Shack-Hartmann wavefront sensor. 

In the analytical investigation presented here we only considered one singularity. In the 
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numerical simulation we found that when different oppositely charged singularities are in 
close proximity to each other, their respective topological charge distributions, as produced 
by the circulation computations, would overlap, causing partial cancelation and a subsequent 
reduction in their peaks. This makes detection of these singularities more difficult. Further 
investigation is needed to understand the effect of multiple singularities located near each 
other on the detection process. 

The phase functions of scintillated optical beams are in general noisy. The noise becomes 
larger as the scintillation increases. The analysis that is provided here does not specifically 
include such noise. It is expected that the ability to detect singularities would deteriorate as 
the noise is increased. This is an important aspect that still needs to be investigated. 



A. Circulation function for a general noncanonical singularity at an arbitrary 
location 

Here we provide the general expression obtained when we compute the circulation, Eq. flTT]) . 
of the sampled average slope G, Eq. ffTOj) . obtained as output from the SHWS when the 
input is the phase function of a noncanonical singularity, located at fi = xq/w and u = hq/w 
with morphology angles a and /?, 
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where Ap = 1 + A and = 1 - A, with A, 5 and C as defined in Eqs. f fT8ll20|) . 
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List of Figure Captions 

Fig. [T] One-dimensional representation of a Sliack-Hartmann wavefront sensor. 

Fig. [2] An array of subapertures (small squares) within the system aperture of the SHWS. 
The average phase slope values are associated with the dots inside the small squares. 

Fig. [3] One-dimensional representation of one lenslet in the SHWS, showing the shift of 
the focal point due to the average tilt of the incident wavefront. 

Fig. m Circulation D"^'"^ over four subapertures with a singularity located either at the 
center (assumed to be the origin), denoted by A or at some arbitrary location {xo,yo) 
denoted by B. The four subapertures are represented by the four squares. The dot at the 
center of each subaperture is the position with which the average phase slope value G of 
that subaperture is associated. The arrows represent the components of G. The dashed-lines 
represent the contour used for calculation of the circulation. 

Fig. [5] Circulation D for a canonical singularity. A topview of D is shown in (a) as a 
function of the relative position of the singularity inside the four subaperture area, shown 
in Fig. m for —2 < /i < 2 and —2 < u < 2. One-dimensional functions of D are plotted as 
functions of r = \fiiF^Pu^ in (b) along the 'diagonal line' and 'middle hne,' respectively, as 
indicated in (a). 

Fig. [6] Peak value of the circulation Z) as a function of the morphology angles < a < vr 
and < /3 < 2t\. The jump at a = 7r/2 is due to the change of the topological charge of the 
singularity. 

Fig. [7] Circulation D for a noncanonical singularity, with a = 7r/4 and /? = tt. A topview 
of D is shown in (a) as a function of the relative position of the singularity inside the four 
subaperture area, shown in Fig. HI for —2 < fi < 2 and —2 < u < 2. One- dimensional 
functions of D are plotted as functions of r = y/ fi'^ + u'^ in (b) along the 'diagonal line', 
'/X-line' and 'i/-line,' respectively, as indicated in (a). 

Fig. [8] Circulation D for a noncanonical singularity, with a = 47r/9 and (3 = tx /2. A topview 
of D is shown in (a) as a function of the relative position of the singularity inside the four 
subaperture area, shown in Fig. HI for —2 < fi < 2 and —2 < u < 2. One-dimensional 
functions of D are plotted as functions of r = in (b) along the 'diagonal line I' 

(perpendicular to the orientation of the singularity), 'diagonal line IF (along the orientation 
of the singularity) and 'middle line,' respectively, as indicated in (a). 

Fig. [9] Numerical simulation results for a Gaussian beam that propagated over a distance of 
100 km through a turbulent atmosphere. The resulting phase of the beam inside the system 
aperture is shown in (a). There are two pairs of oppositely charged phase singularities. The 
pairs are, respectively, located at the lower left and the upper right of the system aperture. 
The circulation D, numerically calculated from the output of the Shack-Hartmann wavefront 
sensor, is shown in (b). 
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Fig. 1. One-dimensional representation of a Shack-Hartmann wavefront sensor. 
ChenFigl.eps. 
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System 
aperture 



iubaperture 
(m,n) 



Fig. 2. An array of subapertures (small squares) within the system aperture of 
the SHWS. The average phase slope values are associated with the dots inside 
the small squares. ChenFig2.eps. 
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Fig. 3. One-dimensional representation of one lenslet in the SHWS, showing 
the shift of the focal point due to the average tilt of the incident wavefront. 
ChenFig3.eps. 
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Fig. 4. Circulation _D™'" over four subapertures with a singularity located 
either at the center (assumed to be the origin), denoted by A or at some 
arbitrary location (xq, yo) denoted by B. The four subapertures are represented 
by the four squares. The dot at the center of each subaperture is the position 
with which the average phase slope value G of that subaperture is associated. 
The arrows represent the components of G. The dashed-lines represent the 
contour used for the calculation of the circulation. ChenFig4.eps. 
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(a) 




Fig. 5. (Color online) Circulation D for a canonical singularity. A topview of 
D is shown in (a) as a function of the relative position of the singularity inside 
the four subaperture area, shown in Fig. HI for —2 < fi < 2 and —2 < u < 2. 
One-dimensional functions of D are plotted as functions of r = ^/J^^~+l^ in 
(b) along the 'diagonal line' and 'middle line,' respectively, as indicated in (a). 
ChenFig5.eps. 
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Fig. 6. Peak value of the circulation D as a function of the morphology angles 
< a < vr and < jS < 27r. The jump at a = 7r/2 is due to the change of the 
topological charge of the singularity. ChenFig6.eps. 
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Fig. 7. (Color online) Circulation D for a noncanonical singularity, with a = 
IT /A and /5 = vr. A topview of D is shown in (a) as a function of the relative 
position of the singularity inside the four subaperture area, shown in Fig. IU for 
—2 < fi < 2 and —2 < u < 2. One-dimensional functions of D are plotted as 
functions of r = v^/i^ + z/^ in (b) along the 'diagonal line', '/x-line' and '//-line,' 
respectively, as indicated in (a). ChenFig7.eps. 
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Fig. 8. (Color online) Circulation D for a noncanonical singularity, with a = 
An/d and (3 = 7r/2. A topview of D is shown in (a) as a function of the 
relative position of the singularity inside the four subaperture area, shown in 
Fig. m for —2 < /i < 2 and —2 < u < 2. One-dimensional functions of D 
are plotted as functions of r = ^/J]P~+l^ in (b) along the 'diagonal line I' 
(perpendicular to the orientation of the singularity), 'diagonal line IF (along 
the orientation of the singularity) and 'middle line,' respectively, as indicated 
in (a). ChenFig8.eps. 
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Fig. 9. Numerical simulation results for a Gaussian beam that propagated over 
a distance of 100 km through a turbulent atmosphere. The resulting phase of 
the beam inside the system aperture is shown in (a). There are two pairs of 
oppositely charged phase singularities. The pairs are, respectively, located at 
the lower left and the upper right of the system aperture. The circulation 
D, numerically calculated from the output of the Shack-Hartmann wavefront 
sensor, is shown in (b). ChenFig9.eps. 
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